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Abstract 

We present and test a new, special-relativistic formulation of Smoothed Particle 
Hydrodynamics (SPH). Our approach benefits from several improvements with re- 
spect to earlier relativistic SPH formulations. It is self-consistently derived from the 
Lagrangian of an ideal fluid and accounts for the terms that stem from non-constant 
smoothing lengths, usually called "grad-h terms". In our approach, we evolve the 
canonical momentum and the canonical energy per baryon and thus circumvent 
some of the problems that have plagued earlier formulations of relativistic SPH. 
We further use a much improved artificial viscosity prescription which uses the ex- 
treme local eigenvalues of the Euler equations and triggers selectively on a) shocks 
and b) velocity noise. The shock trigger accurately monitors the relative density 
slope and uses it to fine-tune the amount of artificial viscosity that is applied. This 
procedure substantially sharpens shock fronts while still avoiding post-shock noise. 
If not triggered, the viscosity parameter of each particle decays to zero. None of 
these viscosity triggers is specific to special relativity, both could also be applied in 
Newtonian SPH. 

The performance of the new scheme is explored in a large variety of benchmark 
tests where it delivers excellent results. Generally, the grad-h terms deliver minor, 
though worthwhile, improvements. As expected for a Lagrangian method, it per- 
forms close to perfect in supersonic advection tests, but also in strong relativistic 
shocks, usually considered a particular challenge for SPH, the method yields con- 
vincing results. For example, due to its perfect conservation properties, it is able 
to handle Lorentz-factors as large as 7 = 50 000 in the so-called wall shock test. 
Moreover, we find convincing results in a rarely shown, but challenging test that 
involves so-called relativistic simple waves and also in multi-dimensional shock tube 
tests. 

Key words: computational fluid dynamics, shocks, special relativity, smoothed 
particle hydrodynamics 
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1 Introduction 



Special-relativistic hydrodynamics has important applications in the fields of 
heavy ion collisions and in astrophysics. Astrophysical examples that involve 
highly relativistic motion include jets from Active Galactic Nuclei [6], pulsar 
winds [18] or gamma-ray bursts [41] and often involve Lorentz factors substan- 
tially in excess of 7 = 10. Analytical solutions are only known for a small set of 
specific problems, for most relevant cases numerical approaches are required. 
A robust special-relativistic scheme can be directly applied to problems of 
the above type and, in addition, it is a core ingredient for general-relativistic 
codes, tools that open up the possibility to tackle a whole new class of inter- 
esting astrophysical phenomena. 

In recent years, grid-based methods for special-relativistic hydrodynamics have 
seen a huge leap forward, and by now many highly accurate Eulerian schemes 
exist, see [30] for a review. Most of these schemes are fine-tuned to solve ID 
relativistic shock problems without oscillations and with sharp discontinu- 
ities. For many astrophysical problems, however, additional capabilities such 
as the accurate advection of smooth flow features are required. In particular, 
for some of the future applications that we have in mind a purely Lagrangian 
method possesses distinct advantages and this is why we focus here on a re- 
fined special-relativistic formulation of the Smoothed Particle Hydrodynamics 
(SPH) method. 

SPH is a Lagrangian, purely mesh-free particle method. Since its first for- 
mulations in an astrophysical context [24,19] SPH has undergone a slew of 
technical improvements and it has found its way into many other branches 
of computationally oriented areas of science. For reviews of the method see 
[5,32,35,45]. Being entirely Lagrangian, the method has obvious advantages in 
advection problems, on the other hand, strong shocks have traditionally posed 
serious challenges. Our aim is to devise a special-relativistic SPH formulation 
that, at least at decent resolution, yields accurate shock-results while keeping 
the other benefits of a Lagrangian scheme. This formulation is intended to 
become the condensation nucleus for a future, fixed-metric implementation of 
general relativistic SPH, the corresponding equations can also be derived from 
a variational principle [36,46]. 

The paper is organized as follows. In Section 2 we derive a SPH formula- 
tion consistently from the Lagrangian of an ideal fluid and the first law of 
thermodynamics and we present the details of our treatment of artificial dis- 
sipation. In Section 3 we investigate the performance of the new equation set 
in a number of special-relativistic benchmark tests. They are complemented 
by the tests shown in [47]. The main results will be summarized in Section 4. 
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2 Special-relativistic SPH with grad-h terms 



In the SPH discretization process, derivatives are expressed as sums over par- 
ticle properties, weighted with the gradient of a smoothing kernel whose width 
is determined by the so-called smoothing length. If symmetrized appropriately, 
the SPH-discretized fluid equations conserve mass, energy, linear and angular 
momentum by construction. In early SPH formulations the derivatives of the 
kernel functions with respect to the smoothing length were assumed to van- 
ish. In practice, however, the smoothing lengths were still evolved to ensure a 
local adaptivity and this inconsistency lead to a violation of the conservation 
properties. How severe this violation is in practice, depends on the consid- 
ered problem [52,43,49]. This deficiency was first addressed by [39] and, more 
recently, by [52] and [34] who derived the SPH equations from a discretized 
fluid Lagrangian. The latter two approaches yielded correction factors for the 
kernel gradients, the so-called "grad-h terms" . 

In contrast to earlier relativistic SPH formulations [21,25,26,22,8,51] we de- 
rive our equation set from a variational principle, similar to [36], but we also 
account for the kernel derivatives with respect to the smoothing length. For 
the flat-space metric tensor, r]^ u , we use the signature (-, +,+,+), Latin indices 
run over (1,2,3), Greek ones run from to 3 with the zero component being 
time. We apply the Einstein sum convention and use c = 1 unless otherwise 
noted. With these conventions the four- velocity, £/ M = dx^/dr is normalized 
to UyP* = -1. 

In labeling the SPH particles we adhere to the following convention: the par- 
ticle of interest is always labeled a and neighbor particles, e.g. those in the 
sum of Eq. (7), are usually denoted by b. If some expression applies to both 
particles a and b of a previous expression, we use the index k. 

2.1 The Lagrangian 

The Lagrangian of a perfect fluid can be written as [17] 



denotes the energy momentum tensor, n is the baryon number density in the 
local fluid rest frame, u is the thermal energy per baryon, s the specific entropy 
and P the pressure. All these quantities are measured in the local rest frame 




(1) 



where 



T» u = (n[l + u(n, s)] + P)CW + Pt\ 



(2) 
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of each fluid element, energies are measured in units of the baryon rest mass 
energy 1 , m c 2 . By using the normalization of the four- velocity, the Lagrangian 
simplifies to 

L pf , sr = - J n(l + u) dV. (3) 



In the general fluid element moves with respect to the frame in which 

the computations are performed ("computing frame", CF). Therefore, the 
baryon number density in the CF, N, is related to the local fluid rest frame 
via a Lorentz contraction 

N = 7n, (4) 



where 7 is the Lorentz factor of the fluid element as measured in the CF. The 
simulation volume in the CF can be subdivided into volume elements such 
that each element b contains u h baryons 

AV4=£ (5) 



These volume elements are used in the SPH discretization process to approx- 
imate a quantity / given at a set of discrete points ("particles") labeled by 
b: 

/(0=£/6^(|r-r 6 U), (6) 
b 7V & 



where our notation does not distinguish between the approximated values 
(the / on the LHS) and the values at the particle positions (/& on the RHS). 
The quantity h is the smoothing length that characterizes the width of the 
smoothing kernel W. The discretization prescription, Eq. (6), yields for the 
baryon number density in the computing frame: 

N(r) = J2"bW(\r-r b \,h). (7) 

b 



This equation takes over the role of the usual density summation of non- 
relativistic SPH, p{r) = J2b m bW(\r — f b \, h). Since we keep the baryon num- 
bers associated with each SPH particle, fix, there is no need to evolve a 



The appropriate mass tuq obviously depends on the ratio of neutrons to protons, 
i.e. on the nuclear composition of the considered fluid. 
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continuity equation and baryon number is conserved by construction. If de- 
sired, the continuity equation can be solved though, see e.g. [8]. The discretized 
fluid Lagrangian reads 

^SPH.sr = TrM 1 + u{n b ,s b )}, (8) 

h iV £> 



or, by use of Eq. (4) 

^SPH.sr = - — t 1 + U ( n b> S b)} ■ ( 9 ) 
b 7& 



2.2 The momentum equation 



The momentum evolution of a particle a follows from the Euler-Lagrange 
equations 

d (9L S PH,sr _ <9L S PH,sr _ g qqx 

dt dv a dv a 



In calculating the baryon number density of particle, b, we use 6's own smooth- 
ing length 

N b = Y,VkW(\r b -r k \,h b ) (H) 
k 



and adapt the smoothing length according to 

l/D 

J 



where 77 is a suitably chosen numerical constant, usually chosen around 1.5, and 
D is the number of spatial dimensions. Hence, similar to the non-relativistic 
case [52,34], the density and the smoothing length mutually depend on each 
other and an iteration is required to obtain a self-consistent solution for both. 
The density gradient with respect to particle position a is given by 



V a N b 



f dW bk (h b ) dr bk dW bk {h b ) dh b dN b 

l^V k \ 5 — + 



dr 



bk 



df a 



dh b dN b df a 



= ttJ2 Vk^ bW bk (h b )(8 ba - 6^), 

"6 k 



(13) 
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where the "grad-h" correction factor 
o - i dh b 9W bk (h b ) 

^= l - m ^^KT (14) 



was introduced. Similarly, the time derivative becomes 



dN a ^ fdW ab (h a ) dr ab dW ab (h a ) dh a dN a 



fa \- ( 



dt j y dr ab dt dh a dN a dt 

= ^J2 ly bVabVaW ab (h a ). (15) 
"a b 

The canonical momentum is given by 



_ <9L S p H , S r d fl + u(n b ,s b )\ , P a\ n ^ 

Pa = K=i = - U b^T = ValaVa U + U a H (16) 

dv a V dv * V 7b 7 V Wo/ 

where we have used the first law of thermodynamics, 

®).-3 



and the relation between the baryon number densities in the different frames, 
Eq. (4). The last term in brackets on the RHS of Eq. (16) is the enthalpy per 
baryon. As numerical variable, we evolve the relativistic canonical momentum 
per baryon, 

S a = la V a (l + U a + ^ ). (18) 



To find its evolution equation dL/df a needs to be calculated. By once more 
using the chain rule, the first law of thermodynamics, Eq. (17), Eq. (4), 
Eq. (13) and V b W ba = —V a W ab , which follows from the choice of a radial 
kernel, W(r) = W(\r\), one finds 



dL 



SPH,sr 



= -E 



v b du 



lb dr, 



7i : n b 



= -v a £ v b (^-V a W ab (h a ) + -^V a W ab {h b ) 

b V*a"a 



P. 



Ngn b 



(19) 
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so that our special-relativistic momentum equation reads 

= - £ (t§t V.W + T^-V.W^^)) . (20) 



£5 T/ie energy equation 



We use the canonical energy to identify a suitable energy variable. We find 

1 + u a \ 



a 



'a - L = "a ( 
a \ 



la 



(21) 



which can be transformed into 

Pa\_Pa 

nj N n 



7o 1 + «o + 



(22) 



As numerical energy variable we choose the canonical energy per baryon 



e a = la 1 + «o + 



n a J N a 



V a - S a + 



1 +Ug 

la 



(23) 



By using Eq. (4) once more one finds 



Ol f l+ Ug \ _ Pg OlNg -> dVg 

dt\ 7a J N% dt a ' dt 



(24) 



and therefore 



de a _ d f - 1 + Mq 



_ dSg P^dNg 

Va ' dt N? dt ' 



(25) 



By inserting Eqs. (15) and (20) into (25), the energy equation becomes 

V a W ab (h b )) , (26) 



de n ^ / P a V b ^ tx t n. \ , P bV a 



dt 



N b 2 Q b 



similar to the non-relativistic case when the "thermokinetic energy" \v 2 + u 
is evolved, e.g. [35]. Alternatively, one could evolve the specific entropy, for a 
discussion see [52,38]. 
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2.4 Artificial dissipation 



Our main aim is an accurate description of an ideal fluid without dissipation. 
We do require however local artificial dissipation to produce entropy at shocks 
to ensure the proper jump conditions, very similar to what nature does on 
scales well below the numerical resolution scale. In that sense one can think 
of both artificial viscosity and Riemann solvers as a subgrid model for physi- 
cal viscosity that would act on an unresolvable scale. Riemann solvers can be 
successfully used in SPH [20,7], but often one prefers a shock treatment via 
artificial viscosity which does not require the restriction to an ideal gas and 
that avoids the explicit solution of the Riemann problem. Guided by the suc- 
cesses of relativistic Riemann solvers [28], Monaghan has constructed a new 
form of artificial viscosity terms [33]. We start from this form of dissipative 
terms, but augment it by a new form of signal velocity and two triggers that 
indicated when to apply it. 



2.4-1 The form of the dissipative terms 

The dissipative terms used in this work are similar to the suggestion of Chow 
and Monaghan [8] 

= -Y^^abVaWa^ with U ab = - (S* a - S* b ) ■ e ab (27) 




and 




-e* b )e ab . (28) 



Here K is a numerical constant of order unity, v s \ g an appropriately chosen 
signal velocity, see below, N ab = (N a + N b )/2, and 

tab = (29) 

\r a - n\ 



is the unit vector pointing from particle b to particle a. For the symmetrized 
kernel gradient we use 

1 



V a W ab = - [V a W ab (h a ) + V a W ab (h b )] . (30) 
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Note that in [8] V a W a b(h a b) was used instead of our choice V a W a b, in practice 
we find the differences between the two symmetrizations negligible. The stars 
at the variables in Eqs. (27) and (28) indicate that in Eqs. (18) and (23) the 
projected Lorentz factors 



ll = i (31) 

yJl-(Vk- Cab) 2 



are used instead of the normal Lorentz factor. This projection onto the line 
connecting particle a and b has been chosen to guarantee that the viscous 
dissipation is positive definite [8]. 



2.4-2 Signal velocity 

The signal velocity that enters the artificial dissipation terms is an estimate for 
the speed of approach of a signal sent from particle a to particle b. The idea is 
to have a physically sound estimate that does not require much computational 
effort. In [33,8] the numerical solution of test problems was found to be rather 
insensitive to the exact form of v sig . For our formulation, we use 

^sig,ab = max(« a , a b ), (32) 



where 

a±=max(0,±A±) (33) 



with A^ being the extreme local eigenvalues of the Euler equations, see e.g. 
[30], 



\$ = — ^- ^ " — (34) 

1 - v 2 d 



and c s fc being the relativistic sound velocity of particle k. In 1 D, this simply 
reduces to the usual velocity addition law, A^ = (vk ± c Sj fc)/(l ± VkC Sy k)- The 
results are not particularly sensitive to the exact form of the signal velocity, 
but in experiments we find that Eq. (32) yields somewhat crisper shock fronts 
and less smeared contact discontinuities (for the same value of K) than the 
suggestions of [8]. 
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2.4-3 Controlling the amount of dissipation 



To ensure artificial viscosity does not influence the flow away from shocks, 
we make the involved viscosity parameters time dependent, a strategy sug- 
gested by [37] and subsequently successfully applied in several approaches, 
e.g. [48,11,49]. The art consists in finding triggers that indicate in which spe- 
cific portion of the flow, or more accurately, at which SPH particle artificial 
dissipation is needed. It is needed both at the shock front itself and in addition 
possibly in the post-shock region to damp velocity noise. What complicates 
things further is that both effects may need different amounts of viscosity, so 
that choosing a large viscosity that is able to resolve a strong shock maybe 
more than what is needed to damp post-shock noise. We therefore aim for 
two independent triggers: one that indicates a shock and another that triggers 
on noise in the velocity field. We follow here the recent suggestion of [9] to 
jump immediately to the desired value of the viscosity parameter rather than 
including the triggers in a source term [37,48] that leads to a continuous rise 
of the viscosity parameter (which in some situations was found to be too slow 
[9])- 

For a particle a we determine a "desired" value of the viscosity parameter, 
see Eqs. (27) and (28), due to a possible shock, if a)S hock, and one due to the 
possible presence of velocity noise, i^noise- The desired value is then 

Ka,des = max (-Ka, shock; Ka, noise)- (35) 



If K ades > K a (t), we instantaneously set K a = K a(ies , otherwise K a (t) smoothly 
decays according to 

dK a K a (t) - K min 

dt T„. 



where 



Xh a 
rmn b (v sig>ab ) 



is the decay time. We use x = 20 for a relatively slow decay. In the earlier 
approaches [37,48] the parameter K min was held at a small but finite value to 
keep particles well-ordered. Since our scheme reacts immediately on noise, we 
can set K min = 0. Our experiments did not show any difference between a zero 
and a small non-zero value. 

Similar to [4,9], this scheme could be further augmented by applying "limiters" 
that suppress viscosity in cases where it is not desirable, but this is beyond 
the scope of the current paper. 
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Shock trigger 

We use the local temporal change of the velocity divergence, d(V ■ v)/dt, to 
identify the emergence of shocks. This is different from earlier approaches 
[37,48] which in their source terms trigger on V • v rather than on its temporal 
change. As already noted in the original paper [37] such a scheme would also 
spuriously trigger on a constant slow compression with V • v— const. Numer- 
ically, we calculate the divergence via "linearly exact derivatives", see below, 
and the temporal change by comparison with the last time step. Note that 
d( V • v)/dt is not actually used to determine the desired dissipation parame- 
ter value, it merely serves as indicator where to take action. This is different 
from [9] who used it to determine a desired viscosity parameter 2 . Instead, we 
calculate the desired shock viscosity parameter via the relative change of the 
density N across the kernel, 

Ai.shock — " ' h a , (38) 



where ViV a is again calculated via linearly exact derivatives, see below. Using 
this very local indicator of the density slope at the shock, we calculate the 
desired dissipation parameter at the shock as 



K, 



a, shock 



K n 
0, 



"^-a, shock 



-^ref , shock "h-^a, shock 



for d(V • v)/dt > and V • v < 
else 



(39) 



In experiments, we find very crisp shocks for a reference value A rc f iS h oc k = 0.9 
together with K m3jX = 2, but to be on the safe side we use A re f iS h oc k = 0.8 in 
the following tests. This produced good results in both 1 and 2D. Note also 
that the actually reached peak values of K are usually substantially below 
-R'max, we will show examples in the context of the 2D shocks, see Sec. 3.2.2. 
We now want to address briefly the exact linear gradients. The derivative of 
a quantity A with respect to coordinate / at position r a can be calculated as, 
see e.g. [43,49], 



M lk A k = M A 



(40) 



with 



M lk = 



n -l 



T,Mn-r a ) l V k a W ab 



(41) 



2 In their work, Eq.(14), the viscosity parameter is determined by a ratio of a 
physical and a resolution-dependent quantity. We suspect that at high resolution 
this will produce too small a viscosity parameter. 
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and 



A k = Y,<A b -A a )V k a W ab 

b 



(42) 



The matrix M corrects for effects from the particle distribution, so that 
linear functions are exactly reproduced even for an irregular distribution of 
particles. 

To illustrate the accuracy of different gradient estimates, we perform a simple 
experiment. We distribute SPH particles, once on a hexagonal lattice and once 
via the regularization sweeps described in Sec. 3.2.1, see Figs. 13 (right) and 14 
(middle), and assign them the same density and baryon number and a pressure 
according to their positions P(x) = 1 + (x a — Xi) 2 /(x2 — xi) with x\ = —0.3 
and X2 = —0.1. Subsequently we calculate pressure gradients according to 
("gradient 1") 

(VP) a>1 = £ ^P b V a W ab (h a ), (43) 
b 7V b 



or ("gradient 2") 




or according to the exact linear gradient. The second estimate is just an SPH 
estimate of (VP) a — P a V(l), see Eq. (6), which just subtracts the leading error 
term from Eq. (43). For the case where the particles are located on the hexag- 
onal lattice all estimates yield accurate results and lie on the same straight 
line as they should. For the irregular particle distribution, the first gradient 
estimate produces a substantial scatter (black) around the exact result, see 
left panel Fig. 1. The second gradient estimate (blue) and the exact linear 
gradient result (red) are hard to distinguish by eye, but the latter produces 
errors that are approximately four orders of magnitude smaller (Fig. 1, right 
panel) . 

Noise trigger 

Our aim is to find an additional trigger that indicates "velocity noise" as it 
can appear behind a shock. Such regions are characterized by some particles 
suffering expansion (V-v > 0) while their neighbors feel a compression (V-v < 
0). Therefore, the ratio 

I s s ¥wf < 45 > 
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Fig. 1. Accuracy of gradient estimates for a disordered particle distribution. Left: 
gradient estimates for the different methods. Right: deviations from the exact result. 

can deviate from ±1 in a noisy region since contributions of different sign are 
added up in 5i )0 and therefore such deviations can be used as a noise indicator. 
To remain very local and to avoid an unnecessary smearing of the shock front, 
our summation in Eq. (45) only runs over neighbors within h a (our kernel 
extends to 2h a ). We use 



S 2 ,a 



(46) 



where the quantity 



Si, a 



-S ha if (V • v) a < 
SVa else 



(47) 



should be zero without noise and our noise trigger becomes 



K„ 



A a - 



- A c ■ I A 
^■-rer , noise '-^a ,n 



for S 2 , a > 0.001c s , a //i a 
else 



(4E 



For the noise reference value we use v4 re f noise = 5. The threshold for S 2 , a was 
introduced to avoid triggering on acceptably tiny fluctuations around zero. 
Note that our choice of dissipation parameters is on the "low-viscosity side" 
and sometimes can produce small, but in our opinion acceptable, oscillations. 
This can be cured, of course, by applying more dissipation. 
The quantities K ajShock and -K^shock ar e stored as an accurate indicators of 
whether a particle is in a shock or a noisy region. 
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2.5 Smoothing kernel 



Traditionally, most SPH formulations use the cubic spline (CS) kernel sug- 
gested by [31], 



W CS (q) 



N 
hP 



1 - \q 2 + f q 3 for q < 1 



1(2 -qf 



for 1 < q < 2 
else 



(49) 



with q = j-, D the number of spatial dimensions and iV the normalization 
(2/3 in ID, 10/7-7T in 2D and 1/tt in 3D). It has been shown to yield good 
results over a large variety of test problems. This kernel, however, has the 
known shortcoming that its vanishing derivative at r = allows particles 
to "pair" once they come close enough to each other, say in a shock. Often 
this has no dramatic effect, but it effectively reduces the resolution due to a 
poorer volume sampling by the SPH particles. Recent investigations [44,54] 
find that kernels that are centrally peaked perform better in Kelvin-Helmholtz 
instabilities since they enforce a more regular particle distribution across the 
contact discontinuity. We find very good results in ID with the CS kernel, but 
in 2D tests we also explore the performance of the centrally peaked "Linear 
Quartic" (LIQ) kernel [54] in the form 

F — u for u < x s 

Au 4 + Bu 3 + Cu 2 + Du + E for u < 1 (50) 
else 



WW?) = -tb ' 



with x s = 0.3, A = -1.458, B = 3.790, C = -2.624, D = -0.2915, E = 0.5831 
and F = 0.6500 and u = q/2. The normalization constant iV is 2.962 in 2D 
and 3.947 in 3D. A comparison of both kernels and their derivatives is shown 
in Fig. 2. The results of our 2D shock test 8 favor the LIQ over the CS kernel. 



2.6 Conversion between primitive and numerical variables 



The new numerical variables, N, S and e, obtained from the integration pro- 
cess, need to be converted into the physical quantities 7, u, n and v. We follow 
the strategies of earlier special-relativistic approaches [29,33]: all variables in 
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Fig. 2. Comparison of kernel W (solid) and its derivative, dW/dq, (dashed) for the 
cubic spline (CS) and the linear-quartic (LIQ) kernel. 



the (polytropic) equation of state 
P = (r - l)nu 



(51) 



are expressed as a function of the updated numerical variables and the pressure 
itself. The resulting equation is solved numerically for the new pressure which 
is subsequently used to recover the physical variables. From Eq. (18) and (23) 
one finds 



S 



e + P/N 



(52) 



and thus 

7 = 



^1 -S 2 /(e + P/Ny 



(53) 



Using Eq. (52) and Eq. (4) one can express the specific energy as 

u=- + -1-7 2 -I- 

7 jJM 



(54) 
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With aid of Eqs. (4) and (54) Eq. (51) can be solved for the new pressure P 
that corresponds to the new values of the integrated numerical variables. Once 
P is known, the Lorentz factor can be calculated from Eq. (53), the specific 
energy from Eq. (54) and the velocity from Eq. (52). 



2.7 Time integration 



We use the optimal third-order TVD algorithm to integrate the system of 
ordinary differential equations of the form du/dt = f{u). The solution is 
advanced by one time step At to time t n+1 according to 



m (1) =m" + At/(iT) (55) 

^) = ^T + + ^Atf(uW) (56) 

^ t+1 = 1 u n + 2 u^ + 2 At / (u&) (57) 
3 3 3 

and we choose a simple error estimate to control the time step. Since the 

computing frame number density plays a central role in the discretization 
process we use it to measure the error growth rate: 



/V n+1 -/V™ +1 
|JV "' ls a,RK2 1 



N, 



n+l 



At 



(58) 



where N^ K2 is the density estimate at t n+1 obtained after a second-order 
Runge-Kutta step. The new time step is then chosen as At new = s(e tol /e) 1 / 2 At old . 
For the "safety factor" we use s = 0.8 and for the tolerable error growth rate 
etoi = 5 x 10 -4 . We use this rather conservative time step choice in the tests 
presented below. We find, however, comparable results with a simpler time 
step choice similar to [8], where the time step, At = min a At a , is determined 



by the momentum change according to At a = 0.3yh a /\dS a /dt\. In our exper- 
iments energy and momentum are conserved to about one part in 10 14 . 



2. 8 Reference formulation 



The performance of the new equation set is compared to the formulation of [8] 
which produces the best shock test results of all published SPH formulations 
that we are aware of 3 

3 Note that [8] obtain their density estimate from integration rather than by sum- 
mation. 
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N a = Y J ^W ab {h ab ) (59) 

b 

^ = "E^(^ + ^ + H^cm) V a W ab (h ab ) (60) 
^ = "E^(^ + ^ + JUcm) • V a W ab (h ab ), (61) 
where ft, a & = (h a + h b )/2 and 



n a fe,cM = -^ 5 (S'* - 5 6 *) • e a6 and fi a &,CM = ^T7^« - 4)^ab, (62) 

where they use a fixed value K = 0.5 for their ID tests. For the following 
comparisons we use their first suggestion for the signal velocity 

°a + \ v ab\ °b + \ V ab\ ,1*1 u * - (ro\ 

v s i g ,CM = , , + ■ , + K 6 1, where v a6 = -v afe • e a6 , (63) 

1 -T <-a|^ a fel 1 c br a fel 



their suggested alternative yields very similar results [8]. 



3 Numerical results 



We use equal mass particles in our tests, so that the density information is 
encoded in the particle separation. Since SPH smoothes "discontinuities" over 
a few resolution lengths, we consider it consistent with the spirit of the method 
to start a simulation from initial conditions that are smooth enough to be 
properly resolved by the method. Throughout the test bench, we approximate 
discontinuities in the initial conditions of a function / via Fermi-functions 

/<*) = i4ife + A ' (64) 



where /l and /r are the values to the left and right of the discontinuity located 
at xs and Ax is the characteristic transition length. We use half of the average 
of the left and right interparticle separation for Ax. The issue of smoothed 
initial conditions is more a matter of taste than of technical requirement. A 
comparison between simulations with Ax as described and Ax = shows only 
minor differences, see below. 

Unless otherwise noted, about 3000 particles are shown and a polytropic equa- 
tion of state with an adiabatic exponent specific to each test is used. 
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3.1 Tests in ID 

3.1.1 Test 1: "standard" relativistic shock tube 



i — 1 — i — 1 — i — 1 — i — 1 — i — 1 — i — 1 — i — 1 — i — 1 — r 




Fig. 3. "Standard" relativistic shock tube [30] at t= 0.35: velocity (in units of the 
speed of light; upper left), thermal energy (upper right), pressure (lower left) and 
computing frame number density (lower right). The SPH solution is shown as black 
circles, the exact solution as the red line. 

This mildly relativistic shock tube (7 max ~ 1.4) has become a widespread 
benchmark for relativistic hydrodynamics codes [29,8,50,10,30]. It uses a poly- 
tropic exponent of T — 5/3, vanishing initial velocities everywhere, the left 
state has a pressure Pl = 40/3 and a density N L = 10, while the right state 
is prepared with Pr = 10~ 6 and Nr = 1. 

The SPH result (circles, at t= 0.35) agrees excellently with the exact solution 
(solid line), see Fig. 3. Only the contact discontinuity at x ~ 0.25 is somewhat 
smeared out. A striking difference to earlier SPH results [22,51] is the absence 
of any spike in u and P at the contact discontinuity. This is a result of the 
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form of the dissipative terms, Eqs. (27) and (28). 

To explore the dependence of the results on the various new elements we per- 
form the following low-resolution (450 particles between -0.3 and 0.3) runs: 
i) use the new equation set, ii) the reference equation set of Chow and Mon- 
aghan [8] iii) the new equation set, but O = 1 to explore the importance of 
the "grad-h" -terms and iv) the new equation set, but K = K max = 0.5 (the 
value chosen in [8]) to explore the effect of the time-dependent viscosity pa- 
rameters. With our parameters (K max = 1.2, v4 re f shock = 0.8, A re f noise = 5) the 
dissipation parameter reaches 0.57, so slightly larger than the value chosen in 
[8] to ensure a fair comparison. The results are displayed in Fig. 4. All nu- 
merical parameters have exactly the same values in all cases. Overall we find 
a good agreement between all the equation sets. For a comparison we show 
the density N since the shock-compressed shell is the most difficult structure 
to capture and therefore shows the strongest deviations from the exact so- 
lution. The "grad-h" terms improve the left edge of the rarefaction fan (not 
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Fig. 4. Comparison of different equation sets, zoom into the density peak of a low 
resolution simulation. 

shown in the figure) and sharpen the left edge of the shock-compressed shell 
(black square vs. blue triangle). The time-dependent viscosity parameters are 
substantially reduced behind the shock which allows the density peak level to 
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reach closer to the correct value (black squares vs. green circles). The main 
difference between the suggested and the reference equation set comes from 
the use of the different signal velocity and the time-dependent dissipation pa- 
rameters, the grad-h terms are only a minor, though welcome, improvement. 
We also briefly compare smoothed vs unsmoothed initial conditions, see be- 
ginning of Section 3. To this end we use a low-resolution setup (about 700 
particles) once smoothed using Ax = 0.5(Axi e ft + Ax ri ght) in Eq. (64) and 
once without smoothing, Ax = 0. Here Aa; le f t /A£ ri g ht are the particle spac- 
ings on the left- /right-hand side. The result for the quantity that showed in 
earlier approaches the largest deviations from the exact result, the specific 
energy u, are displayed in Fig. 5. Overall, we find only minor differences. 



3 




x 



Fig. 5. Comparison between smoothed (black) and unsmoothed (red) initial condi- 
tions for the specific energy in test 1. 



3.1.2 Test 2: strong blast 

The following test with initial conditions (N, v, P) h = (1, 0, 1000) and (N, v, P) R = 
(1, 0, 0.01) is a more relativistic variant of a shock tube and was first considered 
by [40] . It poses a severe challenge since relativistic effects compress the post- 
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shock state into a very thin and dense shell. The fluid in the shell moves at a 
velocity of v = 0.96 which corresponds to a Lorentz factor of 7 s heii — 3.6, the 
shock front moves with a velocity of 0.986, i.e. 7 s hock — 6.0. This test has be- 
come a standard benchmark for relativistic schemes [40,14,28,27,29,16,55,8,13,10,3, 
Overall, the numerical solution (shown at t=0.16, 1800 particles) agrees well 
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Fig. 6. Strong blast [30]: velocity (in units of the speed of light; upper left), thermal 
energy (upper right), pressure (lower left) and computing frame number density 
(lower right). The SPH solution is shown as black circles, the exact solution as the 
red line. 

with the exact one, see Fig. 6. In particular, the intermediate states in ve- 
locity and pressure are well-captured. However, this difficult test is a severe 
challenge and the numerical solution is not free of deficiencies. Somewhat large 
smearing in the internal energy occurs at x ~ 0.15, this is a result of using the 
maximum local eigenvalues rather than a proper spectral decomposition. Also, 
the numerical peak density value exceeds the exact one and the shock moves 
at a slightly too large velocity, effects that decrease with increasing numerical 
resolution. In comparison to [8] both these artifacts are substantially reduced, 
but nevertheless still present. 

We again perform a set of test runs (400 particles between -0.5 and 0.5): i) 
use the new suggested equation set, ii) the reference equation set of Chow 
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and Monaghan [8] iii) the new suggested equation set, but ft = 1 to ex- 
plore the importance of the "grad-h" -terms and iv) the new equation set, but 
K = K max = 0.4 (like in test ii) to explore the effect of the time-dependent 
viscosity parameters. The results are displayed in Fig. 7. At the given reso- 
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Fig. 7. Zoom into the density spike of low-resolution simulations with the different 
equation sets: exact solution (solid, red), our new formulation (black squares), the 
new formulation, but = 1 (blue triangles), the new formulation, but K = K m3uX 
(green circles) and the original formulation [8] (dashed). 

lution, the new formulation (black squares) clearly performs best. This is the 
only test where we see a clear improvement of the solution due to the grad- 
h terms (black squares vs. blue triangles). Again, controlling the amount of 
dissipation has a major effect, the schemes with constant dissipation (green 
circles and dashed line) show the least satisfactory performance. 



3.1.3 Test 3: sinus oidally perturbed shock tube 

Following [12], we explore a shock tube test whose initial right density state 
is sinusoidally perturbed: 

(N,v,P) L = (5,0,50) and (N, v, P) R = (2 + 0.3 sin(50x), 0, 5). (65) 



The main goal of this experiment is to test the ability to transport smooth 
structures across discontinuities. For this relatively mild shock only a moderate 
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amount of dissipation is required, we use K max = 0.2. The numerical result at 
t = 0.35 is displayed in Fig. 8 together with the exact solutions of unperturbed 
shock tubes, once with a right hand side density value of 2.3 (dashed red line) 
and once with 1.7 (solid red line). Note the slightly larger shock speed in the 
latter case (71.7 ~ 1.151 vs. 72.3 ~ 1.149). The numerical solution accurately 
reaches the correct levels of the limiting solutions in the shocked shell. 
We perform again test runs to explore the influence of the equation sets and 




Fig. 8. Shock tube test where a relativistic shock propagates into a sinusoidally 
perturbed medium [12]. The SPH solution is shown as black circles, the red lines 
indicate the exact solutions for unperturbed right hand side densities of 1.7 (solid) 
and 2.3 (dashed). 



parameters (400 particles between and 1): i) the new equation set, ii) the 
reference equation set of Chow and Monaghan [8] iii) the new equation set, but 
Q = 1 and iv) the new equation set, but K = K max , where for the constant-in- 
cases (ii and iv) we use the maximum value from run i) (=0.2). The conclusions 
from this test set are similar to the previous tests: the effects of the "grad-h- 
terms" are visible, but small and the effects from the time-dependent viscosity 
parameters are substantially more important. These tendencies are found in 
all of the subsequent numerical experiments. 
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3.1.4 Test 4: Ultra-relativistic wall shock 

In this test cold gas moves relativistically towards a wall. Upon hitting the wall, 
a shock front forms that travels upstream against the inflowing gas leaving 
behind a hot and dense post-shock region with zero velocity. In the ultra- 
relativistic limit (v — > 1) the shock travels at f shock = — (T — 1), i.e. ^shock = 
1/3 for the polytropic exponent T = 4/3 that we use in this test. The post- 
shock values of density, pressure and specific energy are N p = T/(T — l)iVj, 
P p = 'yTNi, Up = 7, where iVj is the initial density. 

We model the reflecting wall as "ghost" particles streaming with opposite 
velocity from the right towards the wall located at x — 1. For this extremely 
strong shock we use K max = 1. For the initial gas velocity we use a value 
as high as v — 0.9999999998 corresponding to a Lorentz factor of 50 000! 
We further use Ni = 1 and a specific energy of Ui = 10~ 5 . The results of 
the numerical calculation (1000 particles) at t=l are shown together with the 
ultra-relativistic limit values in Fig. 9. The agreement between the numerical 
and the exact result is excellent, only the density and the specific energy show 
minor deviations from the exact solution (maximum error in density 1.2 %, 
in specific energy 1.1%) as a result of so-called "wall-heating" [40] near the 
boundary at x=l. 

3.1.5 Test 5: Relativistic advection of a sine wave 

In this test we explore the ability to accurately advect a smooth density pat- 
tern. We choose a sine wave that propagates towards the right through a 
periodic box. Since this test does not involve shocks, we switch off the artifi- 
cial dissipation terms. We use 500 equidistantly placed SPH particles in the 
interval [0,1], enforce periodic boundary conditions and use a polytropic equa- 
tion of state with Y = 4/3. We impose a computing frame number density 
N(x) = N + ^sin(27ra), a constant velocity t>o = 0.997 corresponding to a 
Lorentz factor 7 12.92 and we instantiate a constant pressure correspond- 
ing to Pq = (r — l)noMo, where Uq = N / / y, N = 1 and u = 1. The specific 
energies of the particles are chosen so that each particle has the same pressure 
Pq. 

The advection of this relativistic sine wave is essentially perfect, see Fig. 10: 
the result after as many as 100 box crossings (circles) is indistinguishable 
from the initial setup (red line), neither wave amplitude nor phase have been 
noticeably affected during the evolution. 

3.1.6 Test 6: Relativistic advection of a square wave 

In this test we advect a square wave through the interval [0,1], again using 
periodic boundaries. Due to the involved steep flanks this problem is substan- 
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Fig. 9. Ultra-relativistic wall shock test: an ultra-relativistic, cold fluid with 
(N,v,u) = (1,0.9999999998, 1CT 5 ) moves towards a reflecting wall at x=l. Note 
that the initial velocity corresponds to a Lorentz factor as large as 7 = 50000. The 
SPH result is shown as black circles, the exact result in the ultra-relativistic limit 
as red lines. 

tially more challenging than the previous one. We represent our box-shaped 
number density profile, N(x), numerically as the sum of two Fermi-functions 
transiting from a lower state iVi ow at x\ to a higher state, iVhigh, and back to 
the lower state at X2- 

N( X ) = iV, ow 4- (AW. " Mow) { 1+ex p ( ^) " ' (66) 

where Ax sets the length scale on which the transitions occur. For the numer- 
ical experiments, we use N\ ow = 1.0, iVhigh = 11 and Ax is set to twice the 
particle spacing in the low density region. We use equal mass particles in this 
test, constant pressure throughout the box is instantiated as in the sine wave 
problem above and we impose again a constant initial velocity of v = 0.997. 
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Fig. 10. Advection of a relativistically moving sine wave (Lorentz factor 7 ~ 12.92) 
across a periodic box: the initial condition of the computing frame number density, 
N, is shown as solid, red line, the result after as many as 100 box crossings is 
overplotted as circles. 

The results are displayed in Fig. 11: the black circles show the initial num- 
ber density profile, the results after 5 and 10 box crossings are shown as red 
circles and blue triangles, respectively. Overall the results of this challenging 
test agree very well with the initial state, but after 10 box crossings the flanks 
have slightly softened, and the low density state close to the flanks has been 
slightly increased. 



3.1.7 Test 7: Evolution of a relativistic simple wave 

Here we present results of a challenging test that involves relativistic simple 
waves and is rarely shown in the literature. 

Relativistic simple waves [53,15,23,2,1] are characterized by the spatial and 
temporal constancy of two of the three Riemann invariants for one dimensional 
fluid flows. The three Riemann invariants are the specific entropy, s, and the 
quantities 

J± = ln( 7 + U) ± [ -dp, (67) 
J P 
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Fig. 11. Advection of a relativistically moving square wave (Lorentz factor 
7 ~ 12.92) across a periodic box: the initial condition of the computing frame 
number density, N, is shown as black circles, the result after 5 box crossing times 
(red circles) and 10 box crossing times (blue triangles) are overplotted. 

where c s is the sound speed and U the x-component of the four velocity. Sim- 
ple waves propagating to the right/left are characterized by the constancy of 
s and J_/J + . In this test, a purely compressive initial sine pulse is set up that 
propagates into a static, uniform medium. According to the results of [2], the 
initial velocity pulse will steepen while maintaining its peak velocity until it 
evolves into a relativistic strong shock. From thereon, the wave will dissipate 
and continuously decrease in velocity and in the contrasts of density and in- 
ternal energy (as measured with respect to the initial unperturbed state). 
Our setup and parameter choice closely follows [2]. We use the equation of 
state of a radiation-dominated fluid, p = k(s)p 4 ^ 3 . The initial state consists of 
an unperturbed fluid state denoted by subscript with an overlaid "velocity 
pulse". Like in [2], we choose the sound velocity in the unperturbed state as 
Co = 0.3 and display normalized, dimensionless quantities. For the initial (local 
rest frame) density we choose n = 1 and we use the length of our initial ve- 
locity pulse, l = 100, as characteristic length scale. The normalized quantities 
are denoted with a "-symbol: jx = p/(nolo), where p is the mass coordinate, 
n = n/riQ, u = u/uq, s = s/sq, P = P/Po and t = t/lo. We proceed in the 
following steps: 

• choose a sinusoidal velocity profile with a maximum t> max = 0.7 as a function 
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of fx. The width of the pulse in normalized mass coordinates is /to = 7r. To 
keep the pulse sufficiently far away from the (fixed) boundaries particles 
(at jl — and 13), we choose /t pea k = 3 as mass coordinate of the velocity 
maximum. For the ease of comparison, we plot the results with a constant 
/t-offset so that our initial setup (black curves in Fig. 12), coincides with 
the initial conditions of [2], see their Figure 5. 
• Once v is known, we calculate the sound velocity [23,2] as a function of jx 



1+cqv^ / l±v\ 



1 

2V3 



l+CQV^ / l+v \ 



l 

2VS 



+ 1 



(68) 



and from this the specific energy and rest frame number density 



4(1 -3c s 2 ) u \u J 



• Finally, the computing frame baryon number density, iV = 771, is calculated 
which, in turn, allows to assign positions and baryon numbers. 

We display in Fig. 12 v, 5n/rio = n — 1, Su/uo = u — 1 and Ss/sq = s — 1, where 
the axis limits and output times are chosen as in Fig. 5 of [2]. The specific 
entropy is nowhere used in our SPH formulation, but can be post-processed 
for the chosen equation of state as s = P 3//4 n _1 . 

Generally, we find very good agreement with the results of [2]. The initial sine- 
pulse continuously steepens while keeping its maximum velocity constant (to 
within about 2%) until a relativistic shock forms (at t ~ 1), see Fig. 12 upper 
left. Note in particular that the positions of the shock fronts agree excellently 
with those found in [2]. The wave subsequently dissipates, thereby producing 
a double peak structure in the density contrast (at t = 1.09, upper right) with 
the dip coinciding with the peak in the entropy (lower right) and the steep 
flank at /t = 3.5 in the specific energy (lower left). The entropy shows a small, 
spurious overshoot at the leading edge which is the result dividing P 3 / 4 and n 
at the shock front, but apart from this, the agreement with [2] is very good. 
We only find minor differences. Our simulations are less dissipative, the second 
density peak, for example, still exceeds the leading one at i — 1.09 and is still 
visible at 7.40 while by the same time it has vanished in [2]. Moreover, our 
specific energy at t — 1.09 shows a clear plateau behind the shock while theirs 
is a very smooth peak (which may just be the result of lower resolution), and 
our entropy peak (at ft ~ 3.75) is slightly smaller than theirs (Ss/s = 0.61 
vs. sa 0.67). 
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Fig. 12. Evolution of a relativistic simple wave: starting from a velocity sine pulse 
propagating to the right into a uniform medium, a shock forms that subsequently 
dissipates the wave. Shown are the velocity (upper left), and the contrasts in (rest 
frame) density (upper right), specific energy (lower left) and entropy (lower right). 

3.2 Tests in 2D 



Mult i- dimensional calculations are complicated by the fact that the SPH parti- 
cles have the possibility to pass each other, which can cause imperfect particle 
distributions and numerical noise (mainly in the velocity). And in fact, noise 
is a major concern for multi-dimensional SPH calculations. Our strategy in 
this respect is threefold: a) since they can be a major reason for noise, we take 
particular care to prepare accurate initial conditions, b) we have implemented 
an artificial dissipation scheme that -in addition to shocks- also triggers on 
velocity noise, see Sect. 2.4.3, and iii) we use relatively large values for the 
parameter r] in Eq. (11). 



3.2.1 Initial particle distribution 

We have performed some experiments with the initial particle setup. We ex- 
plored three types of initial particle distributions: i) a perfect equidistant grid, 
Fig. 13, left, ii) a hexagonal lattice, corresponding to the distribution of the 
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centers of close-packed spheres, Fig. 13, right, and iii) a "glass"-like parti- 
cle distribution obtained in a relaxation process. To produce the "glass-like" 
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Fig. 13. Initial particle distributions: grid (left) and hexagonal lattice (right). 

distribution, we proceed in three steps: initially, the particles are distributed 
according to a Sobol quasi-random sequence [42], see Fig. 14, left panel. In a 
second step, we choose an "adjustment time step", At ad j a = 0.5 (h a /c Sta ), and 
perform several sweeps where we use the force law 

b 



which is a very simple discretization of the Euler equation. In each loop, the 
particle positions are updated according to 

Ta =T a +\ f a A4 jja , (71) 



while the maximum value of the force \f a \ is constantly monitored. This is 
done for as long as the maximum force value is decreasing from one sweep to 
the next. Typically, the maximum force is reduced by two orders of magni- 
tude with respect to the initial Sobol sequence and provides a visually more 
regular particle distribution, see middle panel in Fig. 14. Since the particles 
do not move much during the optimization process, we store for each particle 
a (generous) candidate list which is used for all the sweeps. Therefore, this 
procedure is a computationally inexpensive way to drive the particles into a 
regular distribution. To further optimize the distribution, we "relax" this par- 
ticle distribution to an optimal state by applying a large dissipation constant, 
K — 3, to the momentum, but not the energy equation (we want to obtain 
a perfect particle distribution, but not artificially heat the system). During 
this process, particles at the edges and, for shock tubes, near the transition 
between the two states, are kept fix, so that each state relaxes separately. 



30 



In our tests we have found the best results with particles initially placed 
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Fig. 14. Particle distributions: initial, quasi-random Sobol sequence (left), after reg- 
ularization sweeps (middle) and final configuration after relaxation (right). In this 
process the peaked LIQ kernel was used. 

on the hexagonal lattice. The (somewhat artificial) particle distribution on a 
equidistant grid gave good shock tube results at low resolution, but at higher 
resolution produced a lot of velocity noise once the particles leave the grid. 
This effect was more pronounced for the peaked LIQ kernel, which we still 
prefer for the 2D tests since it performed slightly better in the below tests 
than the standard CS kernel. 



3.2.2 Test 8: 2D relativistic shock tube 

This is the 2D version of test 1 shown above, i.e. the initial conditions are 
(N,v,P) L = (10,0,40/3) and (N,v,P) R = (1,0, 10" 6 ). For this test we place 
140 000 particles on a hexagonal lattice between [-0.4,0.4] x [-0.02,0.02] and 
use our standard 2D parameter set: i] = 1.7,K max = 2.,K min = 0., A rcf shock = 
0.8, A re f inoise = 5 together with the LIQ kernel. The results are displayed in 
Fig. 15 with SPH particle properties as black circles and the exact solution 
as red line. Generally, we find a very good agreement with the exact solution, 
only the contact discontinuity is somewhat smeared out, similar to the ID 
case. Some small high-frequency oscillations in the velocity occur behind the 
shock front. They are mainly caused by the particles that under the action of 
the peaked kernel have to change from the low density lattice into a higher 
density lattice. This is illustrated in Fig. 16 for a lower resolution calcula- 
tion (40 000 particles) where we have plotted the re-scaled velocities over the 
particle distribution. The small post-shock irregularities in the velocity are 
clearly related to the transition region behind the shock front at x 0.08. 
They could be further reduced at the expense of applying more dissipation, 
but we consider the current parameter set as a good compromise between low 
dissipation and absence of noticeable post-shock oscillations. 

In the subsequent low-resolution tests (2000 particles), we illustrate the 
performance of the dissipation control scheme. We perform all tests with 
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Fig. 15. 2D relativistic shock tube test (140 000 particles). Values at individual 
particle positions (i.e. no averaging or smoothing has been applied) are shown as 
black circles, the exact solution is indicated by the red line. 

7] = 1.7,K max = 2.,K min = 0., A re f iSh ock = 0.8, A rc f inoise = 5. and the LIQ 
kernel. The results from the full scheme are shown in Fig. 17 as blue squares. 
Note that the parameter K (green circles) remains at zero up to the arrival of 
the shock front where it jumps to values close to unity. In the post-shock region 
it remains around 0.5 (triggered by noise), and vanishes again in the expan- 
sion fan. For comparison, we perform another test with identical parameters 
but with K = 1 = const everywhere (black circles), and one more where we 
switch off the noise trigger (red triangles). Clearly, the new scheme substan- 
tially sharpens shock fronts, essentially without compromising the post-shock 
region. With shock trigger only, the shock is sharp, but substantial post-shock 
oscillations occur. The two independent triggers allow to apply a different 
amount of dissipation to both phenomena. 

We also perform two tests to gauge the influence of the smoothing kernel, in 
the first one, we want to explore to which extent "pairing" of SPH particles oc- 
curs, in the second, we start from an imperfect initial particle distribution and 
explore the influence of the kernel on the resulting velocity noise. The result of 
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Fig. 16. Zoom into the particle distribution around the shock front in a 2D rela- 
tivistic shock tube test. Shown are the particle positions (black dots) and the scaled 
particle velocities, v a /50, as red squares. Initially the particles were distributed on 
a hexagonal lattice. 

the first experiment is shown in Fig. 18 where we zoom into the shock fronts 
of a low resolution test, once with the CS and once with the LIQ kernel. In 
both cases the initial particle distribution is obtained by the above described 
relaxation process. The CS kernel (left) produces many particle pairs (some 
are highlighted by red ellipses) which deteriorates the volume sampling of the 
SPH particles. The LIQ kernel (right), in contrast, produces only temporarily 
very few pairs directly at the shock front, but the post-shock region is again 
sampled very regularly. 

In the second test, we start from an imperfect initial particle distribution (10 
000 SPH particles in 2D) as produced by the regularization sweeps, see above. 
The particles are subsequently assigned the properties that correspond to the 
left and right state and finally evolved, once with the CS and once with the 
LIQ kernel. The results at t — 0.2 are shown in Fig. 19. The imperfect ini- 
tial conditions introduce some scatter in the velocities (left), but the average 
values still agree very well with the exact solution. To assess the performance 
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Fig. 17. Zoom into the shock region of a 2D, low-resolution (2000 particles), rela- 
tivistic shock tube test. The velocities resulting from the new dissipation control 
scheme are shown as blue squares, the used values for the dissipation parameter 
K as green circles. The red triangles denote the results from the same scheme, but 
with the noise trigger switched off. The results obtained with constant dissipation 
parameter K = 1 are shown as black circles. 



of the different kernels in this situation, we plot in Fig. 19, right panel, the 
deviation of the particle velocities, v a , from the exact solution, v ex , 



5v a = \v a -v e 



(72) 



The LIQ kernel (blue triangles) produces noticeably smaller errors than the 
"standard" SPH kernel (red circles). 

Summarizing our experiments with the two kernels, we find that the LIQ ker- 
nel performs slightly better, it produces a substantially more regular particle 
distribution and yields smaller errors for initially noisy particle configurations. 
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Fig. 18. Low-resolution ("glass") comparison of the cubic spline (CS; left) and the 
Linear Quartic (LIQ; right) kernel for a 2D, relativistic shock tube test. The left 
panel shows a substantial fraction of "paired" SPH particles (some of which are 
highlighted with red ellipses). The peaked LIQ kernel (right panel), in contrast, 
only produces a few pairs directly at the shock front, but the post-shock region is 
very regularly sampled again. 




Fig. 19. Low-resolution, 2D shock simulation starting from an imperfect initial par- 
ticle distribution. Left: velocity for the linear quartic kernel (black circles), average 
velocity (blue square) and exact solution (solid red). Right: comparison of the ve- 
locity errors for the same test, once with the CS (red circles) and once with the LIQ 
kernel (blue triangles). 

4 Summary 

We have derived a new set of special-relativistic SPH equations from a varia- 
tional principle. This work differs from [36] in accounting also for the special- 
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relativistic "grad-h" terms, corrections for usually neglected derivatives of the 
smoothing kernels with respect to the resolution lengths. We have used an arti- 
ficial viscosity prescription that is inspired by Riemann solvers [33] . We trigger 
independently on shocks and on velocity noise. Since the dissipation applied 
in shocks is tuned according to the carefully measured density slope, shock 
fronts become substantially sharper than for a constant dissipation parame- 
ter, see Fig. 17. Contrary to earlier approaches [37,48] we do not evolve the 
dissipation parameter continously to the desired value, but instead increase 
it instantaneously, similar to the approach of [9]. If not further triggered, the 
parameter subsequently decays to zero. None of the triggers is specific to spe- 
cial relativity, both could be applied as well in Newtonian SPH. 
We have carefully tested this new approach in a slew of numerical benchmark 
tests. We find that the relativistic grad-h terms increase the accuracy of the 
method, but usually only have a moderate effect. The improvements are gen- 
erally dominated by the new signal velocity and the time-dependent viscosity 
parameters. The new approach yields excellent results in the numerical ex- 
periments. As expected for a purely Lagrangian scheme, it performs close to 
perfect in pure advection problems, even at large Lorentz factors, see test prob- 
lems 5 and 6. What is more, it also yields accurate results even in very strong 
shock tests, which are usually considered a particular challenge for SPH. For 
example, the scheme is able to accurately handle wall shock problems with 
a Lorentz factor as large as 7 = 50000, see test 4. We also perform a rarely 
shown, very challenging test in which a relativistic simple wave steepens into 
a strong shock and subsequently dissipates (test 7). Our numerical results for 
this test are in close agreement with those of the original paper [2]. In a last set 
of tests, we have explored the performance in 2D, relativistic shocks. Again we 
find very good agreement with the exact solutions. In these multi-D tests we 
have also experimented with the peaked linear quartic kernel [54] which yields 
slightly better test results than the most commonly used cubic spline kernel. 
Numerical experiments [47] show that the scheme is second-order accurate for 
smooth flows and first-order accurate if shocks are involved. 
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